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ABSTRACT 

An  automatic,  finite-step  numerical  procedure 
is  described  for  finding  exact  solutions  to  non¬ 
linear  optimal  programming  problems.  The  procedure 
represents  a  unification  and  extension  of  the 
steepes t-descent ,  and  second  variation  techniques. 

The  procedure  requires  the  backward  integration 
of  the  usual  adjoint-vector  differential  equations 
plus  certain  matrix  differential  equations.  These 
integrations  correspond,  in  the  ordinary  calculus, 
to  finding  the  first  and  second  derivatives  of  the 
performance  index  respectively.  The  matrix 
equations  arise  from  an  inhomogeneous  Ricatti  trans¬ 
formation,  which  generates  a  linear  "feedback 
control  law"  that  preserves  the  gradient  histories, 

H  (t)  ,  on  the  next  step  or  permits  changing  them 

by  controlled  amounts,  while  also  changing  terminal 
conditions  by  controlled  amounts.  Thus,  in  a  finite 
number  of  steps,  the  gradient  histories  can  be  made 
identically  zero,  as  required  for  optimality,  and 
the  terminal  conditions  satisfied  exactly.  One  for¬ 
ward  plus  one  backward  sweep,  correspond  to  one  step 
in  the  Newton-Raphson  technique  for  finding  maxima 
and  minima  in  the  ordinary  calculus.* 

As  by-products,  the  procedure  produces  (a)  the 
functions  needed  to  show  that  the  program  is,  or  is 
not,  a  local  maximum  (the  generalized  Jacobi  test) 
and  (b)  the  feedback  gain  programs  for  neighboring 
optimal  paths  to  the  same,  or  a  slightly  different, 
set  of  terminal  conditions. 


the  whole  interval  tQ  _<  t  <_  tj  where  H  is  the 

variational  Hamiltonian  introduced  in  the  next 
section. 

The  final  time,  tj  ,  may  be  given  either 

explicitly  or  implicitly  in  Eqns.  (4),  For 
simplicity  of  presentation,  we  will  first  discuas 
the  case  where  the  final  time  t1  is  given 
exp  1  i  c  i  t  ly . 

CASE  WHERE  FINAL  TIME  IS  GIVEN  EXPLICITLY 

In  the  usual  manner  we  introduce  the  auxiliary 
scalar  functions 

H(X,x,u,t)  -  L(x,u,t)  +  XTf(x,u,t)  (5) 

$ ( v,x ,t)  ■  4>(x,t)  +  vT(Kx,t)  (6) 

where  X(t)  is  an  n-component  vector  of  influence 
functions  and  v  is  a  q-component  constant  vector. 
We  regard  u(t)  as  control  functions  and  v  as 
control  parameters,  and  introduce  a  modified  per¬ 
formance  index  J  where 


J  -  *[v,x(t1),t1]  +  {H [ X ( t )  ,x(t) ,u(t) ,t] 


-  XT(t)x)dt  (7) 

Note  that  when  (2)  and  (4)  are  satisfied,  (7)  is 
identical  with  (1). 


CLASS  OF  PROBLEMS  TREATED 

The  method  is  applicable  to  a  class  of  non¬ 
linear  optimal  programming  problems  where  one  wishes 
to  determine  control  functions  u(t)  in 
tQ  -i  c  -1  so  as  to  minimize  (or  maximize)  a  per¬ 
formance  index  of  the  form 


J  «  4>[x(t1)  ,t  x  J  +  J  L[x(t),u(t),t]dt  (1) 
co 

subject  to  the  constraints 

x  -  f (x(t) ,u(t) , t ]  x(t)  is  an  n-component  (2) 
state  vector 
u(t)  is  an  m-component 
control  vector 

x(t0)  -  x0  (3) 

^[x(tj)  ,tj  ]  "0  ;  i|i  is  a  q-component  (4) 

vector  (q  <  n) 

A  further  restriction  on  the  class  of  problems 

,  2 

treated  in  this  paper  is  that  we  assume  — j  is  a 

3u 

positive-definite  (or  negative-definite)  matrix  over 


This  latter  technique  is  reviewed  briefly  in  the 
Appendix  to  this  paper. 


Necessary  conditions  for  an  extremal  path  are 


(e.g.  see  Ref. 

1) 

•T 

X1  -  -H 

X 

(8) 

0  -  H 

u 

(9) 

xT(tl>  -  *x 

(v-x(ti)’ti!  5  ‘vAch-t, 

(10) 

Suppose  we  arbitrarily  choose  some  control 
functions  u(t)  and  some  control  parameters  v  , 
integrate  Eqns.  (2)  forward  with  initial  conditions 
(3),  and  Eqns.  (8)  backward  with  boundary  conditions 
(10).  In  general,  Eqns.  (4)  and  (9)  will  not  be 
satisfied.  Now,  consider  a  perturbation  around 
this  path: 

6x  ■ 

f  6x  +  f  6u 
x  u 

(11) 

6X  «= 

-H  6x  -  fT6A  -  H  6u 

XX  X  XU 

(12) 

5H  - 
u 

H  6x  +  H  6u  +  fT6X  specified 
ux  uu  u 

(13) 

6x(tQ)  specified 

(14) 

6X(tt)  = 

[<I>  6x  +  ij^dvl  , 

XX  X  t=t 

(15) 

6^  = 

specified 

x  t=tt 

(16) 

We  may  regard  (11)-(16)  as  a  linear. 
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inhomogeneous  two-point  boundary  value  problem  that 
determines  the  functions  6x(t)  ,  6A(t)  ,  6u(t)  , 
and  the  parameters  dv  in  terms  of  specified 
functions  6Hy(t)  and  specified  parameters  <$x(tQ) 

and  6 ^  .  This  is  very  close  to  the  viewpoint  taken 
by  Merriam  (Ref.  2)  and  Kelley,  Kopp ,  and  Moyer 
(Ref.  3). 

To  solve  this  two-point  boundary  value  problem 
we  may  solve  (13)  for  6u(t)  in  terms  of  6x(t)  , 
6X(t)  ,  and  $H  (t)  ,  provided  H  (t)  is  non¬ 
singular:  J  uu 

6u  -  +  H  6x  +  fT6A]  (17) 

uu  u  ux  u 

and,  upon  substituting  (17)  into  (11)  and  (12),  we 
obtain 


where 


6x 

($A 

A 

B 

C 

v 

w 


A  6x  +  B  6A  +  v 
C  6x  -  A^6A  +  w 

f  -  f  H“  Ml 


-f  H" 1  fT 


-H  +  H  H_1H 

XX  XU  uu  ux 


f  H~ 1 6H 
u  uu  u 


-H  H~‘6H 


XU  uu  u 


(13) 

(19) 

(20) 
(21) 
(22) 

(23) 

(24) 


THE  INHOMOGENEOUS  RICATTI  TRANSFORMATION 

In  view  of  Eqns .  (15)  and  (16),  let  us  introduce 
the  following  inhomogeneous  Ricatti  transformation 
(suggested  in  Refs.  4,  5): 

<5A  (t)  -  P(t)6x(t)  +  R(t)dv  +  h(t)  (25) 

6*  -  RT(t)6x(t)  +  Q(t)dv  +  g(t)  (26) 

where  dv  and  are  constant  infinitesimal 

vectors,  P(t)  ,  R(t)  ,  and  Q(t)  are  matrix 
functions,  and  h(t)  and  g(t)  are  vector  functions. 
Now,  differentiate  (17)  and  (18)  with  respect  to 
time : 

6A  -  P  6x  +  P  6x  +  R  dv  +  h  (27) 

0  -  RT6x  +  RTSx  +  Q  dv  +  g  (28) 

Using  (25)  in  (18)  gives 

6x  -  (A+BP) fix  +  BR  dv  +  Bh  +  v  (29) 

Equating  (19)  and  (27)  and  using  (25)  and  (29)  to 
eliminate  6x  and  6A  ,  we  have: 

(C-ATP-PA-PBP-P)6x  -  ( (AT+PB)R+R]dv 


obtain  boundary  conditions  for  P  ,  R  ,  Q  ,  h  , 

?j  ’ 

and 

P  -  -ATP-PA-PBP+C  ; 

P(tI}  '  I*«xi*xx4vT*xx1t-t1 

(32) 

R  -  -(AT+PB)R  ; 

R(V  ■ 

(33) 

Q  -  -RTBR  ; 

Q(t,)  -  0 

(34) 

h  -  -(AT+PB)h-Pv+w; 

h(t,)  -  0 

(35) 

g  -  RT  ( Bh-Pv)  ; 

g(t.)  -  0 

(36) 

Note  that  (32)  is  a  nonlinear  matrix  differential 
equation  (a  matrix  Ricatti  equation),  while  (33)  is 
a  linear  matrix  differential  equation  using  the 
solution  of  (32);  (34)  is  a  matrix  quadrature  using 
the  solution  of  (33);  (35)  is  a  linear  vector 
differential  equation  using  the  solution  of  (32), 
and  (36)  is  a  vector  quadrature  using  the  solution 
of  (35). 

By  integrating  (32) -(36)  backward  along  with 
(8)  from  to  tQ  (a  "backward  sweep")  we 

generate  all  possible  solutions  to  (11)-(13)  that 
satisfy  the  terminal  conditions  (15)-(16)  .  We  may 
think  of  <25)— (26)  as  "boundary  conditions"  at  time 
t  <  tj  that  are  equivalent  to  the  boundary 

conditions  (15) -(16)  at  time  t  ■  t1  .  Thus  the 

boundary  conditions  at  the  terminal  time  are 
"swept"  backward  to  the  initial  time:  a  "forward 
sweep"  then  generates  the  required  particular 
solution  that  also  satisfies  the  initial  conditions 
(14).  This  is  precisely  the  approach  taken  by 
Bryson  and  Frazier  (Ref.  6)  to  solve  the  linear 
smoothing  problem  except  that  the  sweeps  occur  in 
the  opposite  order;  the  "forward  sweep"  is  the 
Kalman-Bucy  filter  which  involves  a  matrix  Ricatti 
equation,  and  the  "backward  sweep"  gives  the 
smoothing  solution  that  satisfies  the  terminal 
conditions. 

After  completing  the  backward  sweep,  the 
required  values  of  dv  in  terms  of  the  desired 
infinitesimal  changes  6Hu(t)  ,  6x(tQ)  ,  and  6tj> 

can  be  obtained  using  (26) : 

dv  -  [Q_1 (6^-g-RT6x) ]  (37) 

c  c0 

Having  these  values  of  dv  ,  we  could,  in  principle, 
substitute  them  into  (29)  and  integrate  these 
equations  forward  with  (25)  and  (17)  to  find  6x(t), 
6A(t)  ,  and  6u(t)  (a  "forward  sweep"). 

Alternatively,  using  (25)  we  could  regard  (17) 
as  a  linear  feedback  law  for  determining  6u(t)  : 

6u(t)  =  -H‘‘(t)([Hux(t)+£^(t)P(t)]6x(t) 


-  [ (AT+PB)h+Pv-u+h]  »  0  <30) 

In  a  similar  fashion,  substitute  (29)  into  (28) : 

(RT+RT(A+BF) ]6x  +  (RTBR+q)dv  +  RT(Bh+v)  +  g  =  0  (31) 

Viewing  (30)  and  (31)  as  identities,  valid  for 
arbitrary  values  of  6x  and  dv  ,  it  follows  that 
the  coefficients  of  6x  and  dv  must  vanish;  this 
yields  differential  equations  for  P  ,  R  ,  Q  ,  h  , 
and  g  .  Also,  if  we  require  that  (30)  and  (31)  be 
equivalent  to  (15)  and  (16)  at  the  terminal  time,  we 


+  f*(t)R(t)dv  -  6Hu(t)  +  fV)g(t))  (38) 

Note  that  dv  in  (38)  may  be  evaluated  at  the 
initial  time  t  **  tQ  as  was  done  in  (37)  ojr  we  may 

evaluate  it  at  several  intermediate  times  in  the 
manner  of  a  sampled-data  feedback  law  0£  we  may 
evaluate  it  continuously  in  the  manner  of  a 
continuous  feedback  law.  If  we  do  evaluate  dv 
continuously,  then  (38)  becomes 
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6u(t)  -  -H-1  {  [H  +fT(P-RQ_1RT)  ]6x  +  ( £TRQ- 1  ] 6 
uu  ux  u  u 

+  [-6Hu+f^(h-RQ-‘g)]}  (38a) 

Now,  the  first  term  in  square  brackets  on  the  right 
hand  side  of  (38a)  is  a  linear  feedback  on  deviations 
6x(t)  from  the  nominal  state  variable  histories  and 
will  keep  6Hu(t)  ■  0  ,  6<i  ■  0  ,  for  5x(tQ)  f  0  . 

The  second  term  in  square  brackets  is  the  forcing 
function  necessary  to  produce  the  desired  changes 
while  holding  fiH^Ct)  ■  0  .  The  third  term  in 

square  brackets  is  the  forcing  function  necessary  to 
produce  the  desired  changes  6H^(t)  while  holding 

6ili  ■  0  ;  it  is  a  linear  functional  of  fiH^Ct)  ,  and 
vanishes  if  6H^(t)  *  0  .  We  could,  therefore, 

integrate  (2)  forward  (a  "forward  sweep") ,  using 
(38)  in 


u(t)  »  “old(t)  +  Su(t) 

(39; 

6x(t)  ”  x(t)  -  xold(t) 

(40) 

The  advantage  of  this  procedure  over  previous 
gradient  procedures  is  that  we  have  separate,  precise 
control  over  the  desired  changes  6Hu(t)  and  . 

By  repeating  this  forward-backward  sweep  several 
times  we  can  bring  H^(t)  and  1^[x(t1),t1]  pre¬ 
cisely  to  zero  while  increasing  the  performance 
index;  the  required  number  of  steps  depends  on  the 
successful  range  of  linearization  of  (11)-(16).  We 
suggest  that  if  N  steps  are  to  be  used,  it  would  be 
reasonable  to  choose 

6H<r)(t)  -  -e(rVr~1)  (t)  (41) 

J*tr)  •  -E(t)*(r'1)(X(t1),tl]  (42) 

where  «  r/N  and  r  is  the  step  number;  in 

this  way,  larger  and  larger  reductions  in  the 
"residuals"  are  taken  each  step  and,  on  the  last 
step,  the  whole  remaining  correction  is  made, 
bringing  Hy(t)  and  precisely  to  zero. 

LOCAL  OPTIMALITY  -  GENERALIZED  JACOBI  TEST 
AND  CONJUGATE  POINTS 


When  we  have  succeeded  in  bringing  Hy(t)  **  0 
and  4>[x(t1),t1]  -  0  ,  we  have  generated  an 

admissible  extremal  path.  For  this  case,  the  feed¬ 
back  law  (38a)  simplifies  to: 


ia  called  a  conjugate  point  to  the  terminal  manifold 
4i[x(t1)  ,t1  ]  «  0  .  An  extremal  path  is  not  an 

optimal  path  if  it  contains  a  conjugate  point  (see 
e.g.  Ref.  4) . 


INTERPRETATION  OF  THE  MATRICES  P.  Q.  AND  R 


Let  us  define  a  return  function.  V(v,u,x,t) 


which  is  the  value  of 
state  x  at  time  t 


J  in  (7)  when  starting  from 
■v.  tj  uaing  the  control 

functions  u(t)  in  (2)  and  the  control  parameters 
v  .  Infinitesimal  variations  away  from  a  given  set 
of  initial  conditions,  6x(t)  ,  and  infinitesimal 
changes  in  the  control  parameters ,  dv  ,  while 
holding  6H  (t)  *  0  ,  will  produce  an  infinitesimal 


change  in  the  return  function,  6V  ,  given  by 


6V 


-[■ 


XT(t),*T[x(t  ),t 


6x(t) 

dv 


+  i[6xT(t),dvT] 


From  (44)  it  is  clear  that 


P(t),R(t) 

6x(t) 

RT(t) ,Q(t) 

dv 

P(t) 


T 

1  (t)  - 

3V 

T  3V 

’  *  '  TZ  • 

3x(t) 

32V 

32V  0 

3x(t) 3x(t) 

»  R  “ 

3v3x(t)  ’  w 

(44) 


(45) 


32V 


From  (26),  or  (45)-(46)  ,  we  can  also  write 


T 

R  (t) 


3ifr 

3x(t) 


,  Q(t) 


a  ip 

3v 


(46) 


and  we  note  these  quantities  are  similar  to  the 
steepest-asccnt  quantities  l^(t)  and  I.,(t) 
of  Bryson  and  Denham  (Ref.  7). 

If  the  path  is  extremal  (Hu(t)  *  0)  ,  and 
satisfies  the  terminal  conditions  (* [x( t j ) , t j ]  -  0) , 

then  V  -  V (x , t)  is  the  optimal  return  function  of 
Hamilton-Jacobi-Bellman  theory  (see  e.g.  Ref.  8). 
Equation  (44)  ,  using  (26)  with  64>  ■  0  ,  g  ■  0  , 
to  eliminate  dv  becomes 


6V  =  1T6x  +  j  6xT(P-RQ_1RT)6x 


(47) 


which  gives  the  infinitesimal  change  in  the  optimal 
return  function  for  infinitesimal  changes  in  the 
initial  conditions  6x(t)  holding  the  final 
conditions  constant  (6^  **  0)  , 


6u(t)  -  -H^*(t)(Hux+f^(P-RQ_1RT)]6x(t)  (43) 

since  6^  -  0  and  6Hu(t)  “  0  implies  that 

v-w“h*g*0  (see  Eqns.  (23)  ,  (24) ,  (35) , 

(36)).  If  the  symmetric  m*m  matrix  H^u(t)  is 

positive  (or  negative)  definite  and  the  symmetric 

n*n  matrix  P-RQ~lRT  is  finite  over  the  semi-open 
interval  <.  t  <  tj  ,  then  (43)  indicates 

6u(t)  -  0  if  6x(tfl)  "  0  and  we  are  assured  that 

we  have  generated  a  path  that  is  at  least  a  local 
optimal  path.  This  is  a  generalized  Jacobi  test; 

if  P-RQ-1RT  becomes  infinite  at  some  point  this 


SUMMARY  FOR  CASE  WHERE  FINAL  TIME 
IS  GIVEN  EXPLICITLY 


(A)  Estimate  the  control  functions  u(t)  and 
integrate  x  ■  f(x,u,t)  forward  with  given  values 
of  x(tQ)  .  Record  the  constants  4»[x(t1)  ,  1 1  ]  , 

and  the  functions  u(t)  ,  x(t)  . 

integrate 

(B)  Estimate  the  control  parameters  v  and  / 

A  =  -fTA  backward  with  A(t.)  =  [  ♦  +vTiJ>  ]  , 

X  1  X  X  L  L  J 

using  u(t)  ,  x(t)  to  evaluate  f  (x(t) ,u(t)  ,t]  . 

T 

Calculate  H  -  X  f  and  its  derivatives  H  ,  H  , 
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as  you  go. 


(54) 


H  ,  H 

UX  XX 


H“^  must  also  be  calculated. 


While  doing  so,  one  can  verify  that  H  is 

uu 

positive  (or  negative)  definite.  If  H  does  not 

uu 

satisfy  the  appropriate  condition,  better  estimates 
of  u(t)  and  v  are  required  in  (A) . 

(C)  Simultaneously  with  (B) ,  integrate  Eqns. 
(32)-(36)  for  P  ,  Q  ,  R  ,  h  ,  and  g  backward, 
using  the  derivatives  of  H  from  (B)  and  ^(t) 

from  (41).  Record  the  forcing  functions 

Hui(t)['{Hu(c)+fu(t)h(t)1  ■  U<t) 
and  the  feedback  gains 

Hui(t)[Hux(t)+fu(t)P(t)1  ■  K(t) 

Huu(t)fu<t)R(t)  ■  L(t) 

(D)  Determine  and  record  the  parameters  dv 
from  (37) ,  i.e. 

dv  -  q-‘(t0)[e*-g(t0)-FT(t0)6x(t0)l 

(E)  Repeat  (A)  using  the  improved  estimates  of 


°(ti)  ■  ‘jHin +L)lt.t 


fill  .  ill  +  111  i  +  ill  • 

Dt  3t  8x  x  3u  U 

u  -  -H'‘(fT*  +H  +fT*  f+H  f) 

UU  U  Xt  Xt  U  XX  ux 

Equations  (17) -(24)  are  still  applicable  but,  in 
view  of  (49)-(51),  the  Inhomogeneous  Rlcattl  trsns- 
formatlon  beginning  at  (25)  must  be  generalized  to 
the  following: 

(t)l  P(t)  ,  R(t)  ,m(t)*l  fix  lh(t)l  (55 


R  (t) ,Q(t)  ,n(t)  dv  +  g(t) 
lmT(t) ,nT(t) ,a(t)  dt  8(t) 


Differentiating  (55)— (57)  with  respect  to  time, 
using  the  fact  that  d^  ,  dfl  ,  dv  ,  and  dt^  are 

constsnts,  we  obtain 

6X  P  ,R  ,ra  6x  P  h  ( 

0  -  RT,Q  ,n  dv  +  RT  6x  +  g  ( 

n  IT  •  _T  • 


U(t)  "  Uold(t)  +  U(t)  “  K(t)[x(t)-xold(t)]  “  L(t)dv 

(F)  Repeat  (B',  (C) ,  and  (D)  using  the 
improved  estimates  v  , 


(G)  Repeat  (E)  and  (F)  until  H  (t)  -  0  , 

1  -  o  .  u 

CASE  WHERE  FINAL  TIME  IS  GIVEN  IMPLICITLY 

If  the  final  time,  t^  ,  is  given  implicitly  in 
(4),  then  it  is  necessary  to  estimate  tj  for  the 

first  forward  sweep,  in  addition  to  u(t)  and  v  . 

A  few  additional  equations  must  be  integrated  on  the 
backward  sweep  in  order  to  determine  the  required 
dt  j  for  the  next  forward  sweep. 

The  development  ia  the  same  as  in  the  previous 
case  through  Eqn.  (10);  at  that  point  an  additional 
necessary  condition  is  required  to  determine  the 
final  time,  namely  the  transversallty  condition 

fl(x(t1),tl]  -  [*t+V+Ljt-t  "  °  (48) 

The  development  is  again  the  same  up  to  Eqns.  (15) 
and  (16)  which  are  replaced  by  the  following: 

«(tj>  ♦„(*,>  .♦Xj.-O,)  6x<ti)  0  <49> 

d+  -  *x(tj)  ,  0  .nUj)  dv  +  0  (50) 


mT(t  )  ■  [a  -H  H  (H  +fTfi  ) 1  , 

1  X  U  UU  UX  U  XX  t"t , 


6X(t2) 

dt^ 

- 

d  n 

Using  (55)  in  (18)  gives 

6x  -  (A+BP) 6x  +  BRdv  +  BMdtj  +  Bh  +  v  (6 

Using  (55)  in  (19),  together  with  (61),  we  can 
eliminate  dX  and  6x  from  (58)-(60),  and  obtain 
three  equations  like  (30)  and  (31)  in  6x  ,  dv  , 
and  dt.  .  These  three  equations  sre  satisfied 


■  choose  P  ,  Q  ,  R  ,  h  ,  and 
and  m  ,  n  ,  a  to  satisfy 

g  to 

m  +  (AT+PB)m  -  0 

(62) 

.  t 

n  ■  -R  Bm 

(63) 

.  j 

a  -  -m  Bm 

(64) 

B  ■  -m^(Bh-fv) 

(65) 

1  Dt  u  uu  u  x  t-t. 


where  the  boundary  conditions  for  m  ,  n  ,  a  are 
given  by  (52)-(54).  Note  (62)  is  the  same  linear 
vector  differential  equation  as  (33)  whereas  (63) 
anc  (64)  are  simply  quadratures. 

If  (62) -(65)  are  Included  in  the  backward 
integration  sweep,  then  it  is  possible  to  solve  for 
both  dv  and  dtj  at  t  -  tQ  ,  using  (56)  and 

(57)  where  desired  values  of  d^  and  dfl  for  the 
next  step  are  introduced.  The  desired  value  of 
6Hu(t)  must  be  used  in  solving  for  h  ,  g  ,  and  8 

from  (35) ,  (36) ,  and  (64. 
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APPENDIX 

THE  NEWTON-RAPHSON  METHOD  AND  ITS  APPLICATIONS 
TO  ORDINARY  CALCULUS  PROBLEMS 

In  this  Appendix  the  Newton-Raphson  method  is 
briefly  stated.  It  will  be  seen  that  the  Newton- 
Raphson  method  applied  to  optimization  problems 
becomes  a  second-order  iterative  scheme  which  can  be 
applied  in  t..e  neighborhood  of  a  non-singular 
optimum  in  order  to  obtain  rapid  convergence. 

The  formulation  of  second  order  steepest-ascent 
methods  may  be  based  upon  a  simple  extension  of  the 
Newton-Raphson  method  used  to  solve  s  set  of 
simultaneous  nonlinear  equations.  Suppose  one 
wishes  to  find  an  n-vector  x  -  (Xj,..^)  such 
that 

P(x)  -  0  P  -  (pi , . .  .pn)  (Al) 

The  Newton-Raphson  method  generates  a  sequence 
(x(°)  ,xU)  , . . .)  by  means  of  an  iterative  relation 
(A2)  . 

xk+1  -  xk  -  [(If)-1  p]  k  (A2) 

x*x 

The  rationale  for  this  is  obtained  by  expanding 
P(x)  in  a  power  series  around  xk  . 

p(xk+dx)  -  p(xk)  +  If  |  dx  +  0(dx2)  (A3) 

Setting  p(xk+dx)  -  0  ,  one  sees  that 

dx  -  (|f)-lp(xk)  +  0(dx2)  (A4) 

b'f  ignoring  second  and  higher  order  terms  on  the 
right  hand  side  of  (A4)  one  obtains  an  estimate  of 
the  error  in  x  within  first  order  accuracy.  Thus 
(A2)  approximates  the  solution  within  a  second  order 

error.  The  method  naturally  assumes  -jjj  to  be 
nonsingular  in  the  region  containing  (xk)  and  the 
solution. 


The  Newton-Raphson  method  may  be  extended  to 
finding  a  local  maximum  of  a  function  of  several 
variables  f(x)  .  If  f  is  continuously 
differentiable,  a  local  maximum  x  is 
characterized  by  being  a  solution  to  the 
following  equations 

f  -  0  i  -  1, . . .  ,n  (A5) 

xi 


Applying  the  Newton-Raphson  method  to  these 
equations,  one  arrives  at  a  second-order  steepest- 
ascent  method  by  merely  identifying  fx  with 
Pi  in  (A2) .  i 

The  method  may  be  readily  extended  to  problems 
with  constraints.  Suppose  the  maximum  of  f  is 
wanted  subject  to  the  added  constraint 

g(x)  -  0  (A6) 


In  place  of  this  problem  one  may  substitute  the 
problem  of  extremlzing  f+Ag  with  respect  to  x 
and  A  as  independent  variables .  This  problem 
has  no  constraints  and  may  be  handled  as  the  first 
case.  An  extremal  is  characterized  by  (A6)  and 


f  +  Ag  *0  (A7) 

x  x 

k  k 

Expanding  (A6)  around  a  nominal  solution  (x  ,A  ) 
one  obtains  the  following  set  of  linear, 
inhomogeneous  equations  to  solve: 


(f  ■ 


.  .  +  (f  +Ag  )dx  +  g  dA 

k  ,k  xx  °xx/  *x 

x  ,  A 


0  "  *  k  +  *x  dx 

X 

Solving  (A8)  yields  corrections  dx  and  dA  , 
the  second  order  steepest-ascent  method  becomes 


(A8) 


and 


k+1 

x 

k+1 


k  -L 

x  +  dx 
Xk  +  dX 


(A9) 


Several  cautions  must  be  exercised.  One  is 
that  dx  must  be  small  in  order  to  guarantee  con¬ 
vergence,  which  implies  that  the  original  error 
should  not  be  too  big.  Secondly,  the  nominal  and 
the  maximum  must  be  non-singular  and  normal.  This 
is  necessary  to  guarantee  the  inversion  of  the  basic 
equations.  The  non-singularity  condition  guarantees 
that  one  can  solve  for  dx  .  The  normality  con¬ 
dition  guarantees  that  one  can  solve  for  dA  . 
Thirdly,  one  should  note  that  the  second-order 
steepest-ascent  method  seeks  out  stationary 
solutions,  regardless  of  whether  they  are  local 
minima,  local  maxima,  or  saddle  points.  In  order 
to  be  sure  that  the  sequence  converges  to  the  de¬ 
sired  extremum,  the  eigenvalues  of  the  second 
derivative  matrix  must  be  checked.  This  can  be  seen 
for  the  problem  without  constraints  by  substituting 
(A2)  with  P  “  f  into  a  power  series  for  f 
around  xk  . 

f(xk+!)  » f(xk)  -  \  fxf;X  +  °(fx>  (aio) 

k+1  k 

In  order  to  guarantee  that  f(x  )  >  f(x  )  ,  it  is 
necessary  to  assume  fxx  <  0  . 
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